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K^ Abstract 

C^ We numerically implement the variational approach for reconstruction 

^5 in the inverse crack and cavity problems developed by one of the authors. 

^^ The method is based on a suitably adapted free- discontinuity problem. Its 

^^ main features are the use of phase-field functions to describe the defects to 

CN be reconstructed and the use of perimeter-like penalizations to regularize 

the ill-posed problem. 

^ The numerical implementation is based on the solution of the corre- 

sponding optimality system by a gradient method. Numerical simulations 
are presented to show the validity of the method. 
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-^ 1 Introduction and setting of the method 

O 

r ■ We consider a homogenous and isotropic conducting body, assumed to be con- 

^^ tained in fl, a bounded, Lipschitz domain of M^, N > 2. We assume that 

lO there exist ili, a Lipschitz domain contained in, and different from, il, and a 

^^ closed set 7 C d^ n dfli such that the interior of 7 is not empty and 7 has 

^— ^ a positive distance from r2\0i. We assume that 7 is known and accessible to 

. . measurements. 

^ In the body there might be present some defects, which we assume to be 

k> perfectly insulating and outside fii. Namely, we model these defects by a closed 

;_i set Kq C il, such that Kq n Qi is empty. We notice that Kg represents the union 

d of the boundaries of these defects and that we denote with Gkq the connected 

component of fl\K() containing J7i, that is the region of fl reachable from 7 

without crossing Kq. 

The defects may have different geometrical properties. For instance, we may 
have, even at the same time, cracks (either interior or surface-breaking), or 
material losses (either interior, that is cavities, or at the boundary). We recall 

*Univcrsity of Graz, Institute of Mathematics, Heinrichstr. 36 8010 Graz, Austria. E-mail: 
wolf gang. ringOuni-graz .at 

tUniversita degli Studi di Trieste, Dipartimento di Matematica c Informatica, via Valerio 
12/1 34127 Trieste, Italy. E-mail: rondi9units.lt 



that a defect Kq is a material loss if Gkq coincides with the interior of its 
closure. 

Let us consider the foUowing experiment. If a current density /o is applied 
on 7, then the electrostatic potential in fl, uq = u{fo, Kq), is the solution to the 
following (normalized) Neumann boundary value problem 

{Au = in n\Ko 

Vu-v = /o on 7 
Vu-v^Q on d{n\KQ)\-i 

The current density is modeled by a function /o G L*(7), for some constant 
s > N — 1, such that / /o = 0. The electrostatic potential uq may then be 
measured on 7. We call such a measurement go = Uq\^ and we observe that 
(?o G L'^il) and J go ~ 0. In this way we obtain an electrostatic boundary 
measurement of voltage, go, and current, /o, type on 7. In mathematical words, 
we measure the Cauchy data {go, fo) of the harmonic function uq on 7. Clearly, 
prescribed the current fo, the voltage go depends on Ko- If Ko is unknown, then 
the measured voltage go may provide information about the unknown defect. 
In fact, the aim of the inverse problem is to reconstruct an unknown defect 
Kq by prescribing one or more current densities /o and measuring the corre- 
sponding value of the potentials on 7. Such a problem arises, for instance, in 
non-destructive evaluation, for the determination of flaws like cracks or cavities 
in conducting bodies by non-invasive methods. We refer to this problem as the 
inverse crack problem, in the general case. Instead, when we a priori know that 
the defect is a material loss, we denote it as the inverse cavity problem,. For re- 
sults on the inverse crack problem and related problems, we refer to the review 
article [51 . Here we simply wish to note that a single measurement (that is per- 
forming the experiment previously described only once) is enough to determine 
uniquely a material loss. In the general crack case, instead, one measurement 
may not be enough, however two suitably chosen measurements (corresponding 
to two suitable prescribed current densities) are enough for unique identification 
of any kind of defects at least in the planar case. 

Let us remark here that if the unknown defect is a priori assumed to be 
interior (that is Ko C fi) and if the whole boundary of fl is accessible, then we 
may simply take 7 — dfl. 

Our approach to this inverse problem is the following. We observe that uq is 
smooth outside Ko, whereas it may, and generally does, jump across Ko- There- 
fore, starting from the Cauchy data, we wish to reconstruct the function uq in 
f2, and in particular its discontinuity set J{uo). We notice that this is not a clas- 
sical Cauchy problem for uq, since Uq is harmonic in fl\Ko with Kq unknown! 
Rather, it looks more like a free-discontinuity problem for uq, since its discon- 
tinuity set J{uo) is unknown and it is actually the aim of our reconstruction. If 
we are able to reconstruct uq and J{uo), then we obtain valuable information 
on Kq, given the fact the J{uq) C Kq- Actually, for the inverse cavity problem, 
J(uo) determines the whole OGko- On the contrary, in the inverse crack prob- 
lem, it may happen that a crack is not visible for a particular measurement, 
that is J{uo) does not detect the whole OGko- In this case, we may change the 
prescribed current density, reconstruct again the electrostatic potential from its 
values on 7, and recover another portion of OGko- The uniqueness results tell 



us how many times and with which kind of prescribed current densities we need 
to repeat this procedure to fully reconstruct the unknown defect. 

The main difficulties in the reconstruction of uq from its Cauchy data are 
the following. First of all, the problem is severely ill-posed, as Cauchy problems 
for elliptic equations are. Second, since the potential uq to be reconstructed is 
a discontinuous function whose discontinuities are unknown (actually they are 
the aim of our reconstruction) , the problem is not even linear. Thus all the main 
difficulties of the original inverse problem are still present in the reconstruction 
of Wq. 

The way to tackle ill-posedness is crucial. In fact, since the boundary data 
are measured, the data which are really available are not the exact Cauchy data 
(<?Oj/o) but some noisy perturbation of them. Namely, the available data we 
assume to know are {g^, fe)- Here fe belongs to L''{dn) and satisfies supp(/e) C 
7 and Jg^ f^ = 0, whereas g^ belongs to ^^(7) and satisfies / g^ = 0. We assume 
that 

(1-2) ||/o-/£||l=(7) < e and Wgo - ge\\m^) < £■ 

where e, < £ < 1/2, denotes the noise level. 

As mentioned, rather than a classical Cauchy problem, we consider such a 
reconstruction as a free-discontinuity problem for the unknown potential uq. 
We follow the variational approach developed in (Hlllj for cracks and material 
losses, respectively. Such a method is based on the following two features. The 
first one is the choice of the regularization. In order to regularize the problem a 
perimeter-like penalization is used. Namely, we penalize the (TV— l)-dimensional 
measure of the unknown defect Kq (actually of the discontinuity set of the 
unknown potential). Second, we model discontinuity sets through phase-field 
functions, thus obtaining a formulation in which a discontinuous function u and 
its discontinuity set J{u) are replaced, respectively, by a smooth function u and 
by a smooth phase-field function v. Such a formulation is amenable to numerical 
implementation. 

In particular, for the crack case, in [8j it has been used a regularization based 
on the so-called Mumford-Shah functional, [7j, and its approximation, in the 
sense of F-convergence, with phase-field functionals due to Ambrosio and Tor- 
torelli, [HI2I. For material losses, [9 , it has been used a more classical perimeter 
penalization and its approximation, again in the sense of F-convergence, with 
phase-field functionals due to Modica and Mortola, [B]. For further details on 
free-discontinuity problems and their approximations we refer for instance to 

mm- 

In this paper we develop the numerics of the approach in [51 [S] . The con- 
vergence analysis done in these papers provides a justification of the numerical 
method, in particular for the material loss case, see [3 Theorem 4.2]. For what 
concerns the crack case, we do not have a precise convergence result for the 
method implemented here. However this is quite simpler from a numerical point 
of view than the one developed in [5] and for which we have convergence results. 
Moreover we believe that this simplification might still lead to good reconstruc- 
tions, see also the discussion in Section 5 of |H|. 

We shall use the following notation. We fix a constant q > 2, depending 
on some regularity properties of the defect Kq to be reconstructed. Namely, we 
assume that there exist a constant q > 2 and a constant C, independent of /o. 



such that 

||Vuo||l9(o) < C'||/o||ls(t-), 



where uq — u{fQ,Ko) solves (1.1). We also set 



0<qi = {q-2)/{2q)<l/2. 

The function ^ ; M — > M is continuous and non-decreasing and such that 
7/1(0) = 0, -0(t) > if t > 0, and ipil) ^ 1. Then, for < e < 1/2, we define 

i;, ^ {I - e^)tlj + e\ 

We introduce a single- well potential V centered at 1, that is a non-negative 
continuous function such that V{t) = if and only if i = 1. 

We shall also need a double- well potential W centered at and 1, that is a 
non- negative continuous function such that W{t) = if and only if f e {0, 1}. 

We assume that the functions W, V and ?/j are C^ and are bounded all 
over M and also their derivatives are bounded and uniformly continuous all over 
M. We shall also assume that W', V and ip' are Holder continuous for some 
exponent a, < a < 1, all over M. Furthermore, the following assumption will 
be made. About ip we require that for any i < we have ^p{t) = ip{0) — and 
i/j(i) = -(/(l) — 1 for any t > 1. In particular, we have that ^^'(0) — '0'(1) — 0. We 
also require that for any t < we have V{t) > V{0). Obviously, we have that, 
for any t < 0, W{t) > W{0), and, for any t > 1, V{t) > V{1) and W{t) > W{1). 

For example, the following choices may be made. For any t e [0, 1] 

^P{t)^-2t^ + M^, W{t)^9t^it-lf, V{t)^{t-lf/A, 

with straightforward extensions beyond [0, 1]. 

We define the space W{n) = {w e W^^^{n) : w = a.e. in fli}. To any 
V £ W{il.) we associate the function v — I — v. We remark that v € VF^'^(ri) 
and V = 1 almost everywhere in J7i. We finally fix positive tuning parameters a, 
b and c, and a noise level e, < e < 1/2. All these constants and the notation 
will be kept fixed throughout the paper. 

Basically, the method is the following. Beginning from the crack case, we 
wish to minimize, with respect to the phase-field variable v g W{fi), with the 
constraint < w < 1, the functional J> : T4^(J7) — ?► K, which is defined as follows. 
For any v S W{n,), recalling that w = 1 — w, we set 



e"^^ ./■y Jo £ Jn Jn 



(1.3) F,{v)^—j \ue^9eV + hj Mv)\yu,\' + - I Viv) + e I \Vv\'. 
Here Ue — Us{v) solves 



div{tps{v)Vu^) =0 in ri 

(1.4) { ljJeiv)'VUe ■ V ^ fe OH OD, 

We notice that the first term is the fidelity term with respect to the measured 
boundary datum, the other three terms are the Ambrosio-Tortorelli functional. 
The link with the prescribed boundary datum and with the presence of cracks is 
through Ue, the solution of the weighted elliptic equation. We observe that the 



single- well potential V forces the phase-field function w = 1 — w to be equal to 1 
except in a small region, which is where the crack should be located. The tuning 
parameters a, b and c allow to put more emphasis on one or the other of the fea- 
tures of the functional. Namely, a controls the match with the Dirichlet datum, b 
the smoothness of the reconstructed potential away from its discontinuities and 
c the penalization on the (TV — l)-dimensional measure of the discontinuities. 
Therefore c may be seen as a regularization parameter. 

For the material loss case, we simply replace the single- well potential V with 
the double-well potential W. Namely, we define Gs '■ W{il) — )• M in an analogous 
way by simply replacing V with W, that is, for any v G W{i^), we set 

(1.5) g,{v) ^^ f\u^-g^\^ + b f Mv)\^Ue\' + - [ W{v) + 8 [ \Wv\\ 



We then minimize, with respect to the phase-field variable ii e M^(f^), with the 
constraint < w < 1, the functional G^. 

We notice that in this case the last two terms are the Modica-Mortola func- 
tional, which penalizes the perimeter of Gug in Q. Here the double- well potential 
W forces the phase-field function v to be either (inside the material loss) or 1 
(outside the material loss), with a quick transition between these two regions. 

Summarizing, we shall minimize the functional J>, when we aim to recon- 
struct defects such as cracks, and the functional Q^, when we aim to reconstruct 
material losses. Namely, we wish to solve numerically the following minimization 
problems (depending on the properties of the unknown defect Kq) 

(i) niin J> on W{Q), with the constraint < w < 1, if Kq contains portions 
of cracks. 

(ii) minQ^ on W{fl), with the constraint < v < I, li Kq is a material loss 
defect. 

Let us notice that, by the direct method, these minimum problems admit a 
solution. 

About the numerical method, in order to find the minimizers, we formulate 
the corresponding optimality system and we use a gradient method, see Sec- 
tion |2] for details. In Section [3] numerical simulations are presented for both the 
single- and double-well approximations. Numerical experiments are performed 
for various types of defects with noise-free and noisy data-sets. 
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2 Optimality system and the gradient method 

We now look towards the numerical implementation of the method. We begin 
by recalling the differentiability properties of the functionals J> and Q;, , which 
have been investigated in [Qj Section 6] . 

We define the following spaces. For any p, 2 < p < +cx), let us call Lp{il) — 
{i e LP{n) : V = a.e. in fii} and Wp{n) = W^'^{n) n Lp{n), with norm 

\\v\\L^{n) = \\v\\Lp{n) and ||i)||wp(n) = l|w||Lp(f2) + W^^vllL^ny To any i e L^{n) 
we as usual associate the function v — 1 — v. li v belongs either to Lp{^) or to 
Wp{n), then V £ LP{U), v ^ I almost everywhere in i7i, and, provided < i; < 1 
almost everywhere in 17, wc also have < i; < 1 almost everywhere in 17. We 
observe that W2{^) = W{VL) as previously defined. 
For any q, q > 2, we define 



W^i'«(17) = |ue W^i'«(17) 
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We observe that, by a generalized Poincare inequality, on W!^'*(r2) the usual 
W^''^{n) norm and the norm ||u||^i,9(q-| — || Vu||^g(j7) are equivalent. Therefore, 
we shall set this second one as the natural norm of VF2''(f7). 
We define He : L2in) -> W]^'^{n) as follows 

'He(i)) = Ue{v) for any v € L2{fl)- 

There exist constants p{e) > 2 and q{e) > 2, depending on s and a, such 
that all the following results hold. 

First, He ■■ L2in) -^ WY'^'-''\n), with bounded image in WY''^^\n), and, 
for any vq G ^2(^)7 such an operator 'He is differentiable in vq with respect 
to the LP{n), with p > p{e), and W^''^^''\n) norms. Let DHeivo) ■ Lp{n) ~> 
W7'' (J7) be the differential in vq. Then for any v in Lp{fl) we have 

where Ue = Ue{vo,v) E Wy^(f7) solves the following problem 

" div(7A,(t;o)VC/,) = div(^[(t;o)W(He(So))) in n 



^^■"^^ \ VeK)VC/e-:/ = on 917. 

Obviously, vq — 1 — v^. We recall that for any vector valued function G G 
L^(17,M^), div(G) defines a functional on VF^'^(17) in the following way 

div(G)[0] = - / G • V0 for any e W^'^{n). 



Therefore, the weak formulation of (2.1 1 is looking for a function [4 G VK^'^(f7) 



such that 



V'e(«o)VC/, • V(p = / i^'^{vo)iW{H,{i>o)) ■ V<^ for any cp e W^^^^^^)^ 
n in 

Here, and analogously in the sequel, the differentiability has to be understood 
in the following sense. For any v in Lp(f7) 

He{io + i)= Heiio) + DHeiioM + R{v) 



where 

n n 1™ iFTi = '^• 

\\v\\LP(n)^0 \\V\\LP{n) 

We conclude that, for any vq E W{ft), J> and Q^ arc difFcrcntiablc in vq with 
respect to the Wp{V,) norm, with p > p{e). Let DJ^s{vo), DGs{vo) : Wp{il.) -^ R 
be the differentials in vq of J> and tj^, respectively. Then, for any v € Wp{fl) 
we have 

(2.2) DT.iioM ^^J{H,{io)-ge)U,ivo,i) + 

b f (27^eK)VHe(5o)-Vt/e({)o,^)-^i(«o)|VHe(i)o)|'i^) + 
Jn 

- / {-V'{vo)v) +2e [ Vio ■ Vu 
and 

(2.3) Dg,{vo)[v] = ^ f{He{io)-ge)Ue{io,v) + 

Jn 

— I {-W'{vo)v) + 2e I \/va ■ Vu. 
£ Jn Jn 

An important remark is the following. If A^ ~ 2, then we may actually 
choose p{e) = 2, and we observe that W2{^) is a Hilbert space, with the scalar 
product /j^ Vwi • Wv2 for any ui, z;2 € PF2(^)- If iV > 2, then it might happen 
that p{£) > 2 and that Wp(g)(il) has not a Hilbert space structure anymore. 
However, since p{e) is finite, Wp(^^-){ft) is still a strictly convex real reflexive 
Banach space. 

In the sequel we shall fix p = p{e), (with p(e) = 2 if A = 2) and we call 
MJ^e the following functional, which is defined on W^'^{il) x VFp(e)(0), 

(2.4) MMu, i)^^J\u- g,\' + J {bMv)\Vu\^ + jV{v) + s\Vv\'), 

for any (u, w) e W^'^{n) x %(e)(f^). 

Such a functional is finite for any {u,v) G Wy^(rj) x Wp(j,){i^). By similar 

reasonings, for any {uo,vo) € VF^'' (fi) x Wp(g)(f2), we have that M.T^ is 
differentiable in (uo,vq) and for any (m, w) S VF^'^(r2) x Wp(g)(r2) we have 

(2.5) DMTe{uo,vo)[{u,i})] = ^ {ua ~ ge)u+ 

{2tl;,{va)Wua-Wu-^'^{vo)\Wuo\^i) + - / {-V\vo)i) + 2e / Vwq-Vw. 

£ Jo Jo 

We observe that J>('5) = AlJ^e('He («),«). Analogously, we define AiG^ sim- 
ply by replacing V with VF. Analogous properties of differentiability hold for 
AiGe as well. 



Let us finally define £J"^ : VF|2(^) ^ Wp(e)(fi) x W^^'^{n) -^ R such that 
for any (w, w, 0) e W%i'2(r2) x Wp(e){^) x T4^i'2(Y}) ^g jjave 

(2.6) CT,{u,v,^)^MT,{u,i)+ f ij,{v)Vu-V(b- f f,^. 

Jn Jdn 

In an analogous way we define CGe replacing MT^ with A4Gs- 

We observe that CTg (and £tjg as well) is differentiable in any (mq, vq,(J)o) G 
m4'«^^^(0) X Wp(e)(rj) X W^^^in). For any (u, «,</.) e M^^^'^d^) ^ Wp(e)(f^) x 
W^^^{n) we have 

(2.7) '' (Mo,-go>o)M = 
aw 

— {uo-9e)u + 2b I V'£(wo)Vuo • Vu+ / ?Ae(wo)V0o • Vm, 
and 



£''' J7 Jf2 Jo 



(2.8) ^(^,o,5o>o)[«]--& /^;(«o)|V^io|'5+- f{-V'ivo)i) + 

2e Wio-Wv- Ve(«o)i5VMo • V(/)o, 
Jo Jo 

and, finally, 

(2.9) — -^(wo,«o,0o)[</']= / ^s(«o)Vuo-V</>- / /,(/.. 

00 Jo Jao 

Then the resulting optimality system is the following. We look for critical 
points, or better minimizers, of J>, or, equivalently, of AiT£{u,v) subject to the 
constraint u = He{v). We use a gradient method, whose algorithm is divided 
into steps. A completely analogous method may be used for finding minimizers 

Of^e- 

Step 0: initialization. 

We initialize the algorithm by putting fc = and choosing an initial guess 
vq G VF(ri) such that < wq < 1 almost everywhere. We observe that taking 
i^o ^ (that is Wo = 1) is not a good choice because this is a critical point of 
the functional J>, thus the gradient method fails in this case. 

Step 1: finding u^. 

We solve 

{div(i/'e(wA:)Vufe) =0 in f7 
V'e(l'fc)VMfc -ly^ fe on on 
I^ Uk = 0, 

that is we look for u^ G T4^y^(i7) such that 

(2.11) / Mvk)yuk -Vcj}- ( fe(j) = for any e W^^^i^). 

Jo Jao 



We noti ce t hat Uk = 'He('Cfc) and Uk actually belongs to W^''^ (Q). By (2.6) 
and by (2.9), we have that for any <}> g W^''^{Q) 



CFe{uk,Vk,4>) ^ MTe{uk,Vk) =Te{vk) and 



dCT, 



{uk,Vk,4>) = 0. 



Step 2: finding (pf,. 

We solve the following boundary value problem 



(2.12) 



div{'ipeivk)^(t>k) = -div(26V£(wfe)Vwfe) in fi 
il}e{vk)V(j)k ■ i^ ^ - — {uk- ge)Xj on dfl 



Here x-y denotes the characteristic function of 7, that is 



(ufc - ge)x-, 



(uk-ge) on 7 







on dfl\"f. 



The weak formulation of (2.12) is looking for (pk G W2' (fi) such that 



(2.13) / MvkW^k-^u = 
Jn 



2a 



2b I 'iJj^{vk)Vuk • Vw - ^^ I [uk- ge)u for any u G W^''^(^). 
Jn ^''^ Jj 

Such a solution (pk exists and is unique. Then CT^iuk, Vk, ipk) = ■MTe{uk, Vk) 



Te{vk) and, by ( p77| ) 
dCF. 



dCT, 



iuk,Vk,(t)k) =0 and — - — {uk,Vk,Pk) = 0- 



Step 3: computing the gradient and updating v^. 

We compute the differential of J> at the point Vk- We observe that if u = Heiv), 
then for any (p g W^'^{il) we have 

T,{v) = MT,{HAv), i) = CTAHe{i),i, P). 
Therefore, since Uk — 'He(wfc), and if we pick (p — pki then 

DF^{vk) = '^ (uk,Vk,pk)- 



We conclude that, by (2.8), we have for any v e Wp(g)(f2) 



(2.14) DT,{ik)[v] = -h / ^;(i;fe)|V«fc|2f} + - / {-V {vk)v)+ 



2e / Vwfc • Vw - / ip'^{vk)v'Vuk ■'Vpk- 
Ja Jn 

Let us now consider the space Wp(g)(fi). We recall that either H/p(£)(f7) = 
W2{fl) (if iV = 2), that is Wp(^^-^{V.) is a Hilbert space, or Wp(£)(51) is a strictly 



convex real reflexive Banach space (if A^ > 2). In either cases, iiW — Wp(j)(ri), 
we fix an operator T : W* — ?► W such that for any w* G W* , we have 

{w*,T{w*)) ^\\w*\\'^ and ||r(ii;*)|| = ||7«*||, 

where (•, •) is the usual duality between W* and W. We may choose T as the 
duality mapping from W* into W** = W. If T4^ is a Hilbert space and we 
also identify W* with W, then T is actually the identity. See, for instance, [101 
Section 42.6]. Let us call T^ the corresponding operator for Wp(£)(ri). 
For a positive constant tk, we then update Vk by setting 

Wfe+i =Vk- tkT^{DFe(vk))- 

We observe the following. If DFe{vk) ~ 0, then {uk,Vk,(j)k) is a critical 
point of CT^ and Vk is a critical point of J> and the algorithm comes to a 
stop. Otherwise, provided tk is small enough, an easy computation shows that 

Teivk+l) < ^e{vk)- 

Step 4: normalization and finding ■Ufe-i-i. 

We normahze v^+i by truncation as follows. We set Wfe+i — (i)fe+i A1)V0. In such 
a way we obtain that ■Cj.+i S WpU\{^) and < "Dfe+i < 1 almost everywhere in 

n. 

Let us note that, by our hypotheses, such a truncation does not increase the 
value of the functional, in fact for any ij € Wp(e)(ri), if w = (i) A 1) V 0, then 

J'e{v)<J'e{v). 

Therefore, we have found that either DF^{vk) — 0, and the algorithm stops, or, 
otherwise, provided tk is small enough, J^^{vk+i) < J^s{vk)- 

Once we have computed Wfc+i, we iterate the algorithm by going back to 
Step 1. 



3 Numerical experiments 

The data for the numerical experiments are generated by solving Laplace equa- 
tion numerically on an domain with certain prescribed defects (cracks or cavi- 
ties). We solve the Neumann problem with given flux on the boundary of the 
computational domain, and read off the corresponding Dirichlet data to get 
a feasible pair of Neumann and Dirichlet boundary data on a discrete set of 
measurement points on the boundary from which the defect has to be recon- 
structed. As input fluxes we choose pairs of plus-shaped current profiles with 
opposite sign located at two different sides of the rectangular computational 
domain. The Laplace equation is solved on a very fine irregular grid using lin- 
ear finite elements. The boundary data are genuinely defined on the unevenly 
distributed nodal points of elements on the boundary and are interpolated onto 
a much courser regular grid of measurement points. When experimenting with 
noisy input data, both boundary values are contaminated by adding Gaussian 
distributed artificial noise to the data, usually with different noise levels for 
/= i^k and5 = u|^. 
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For the numerical implementation of step 1 in the algorithm described in 



the previous section (that is the numerical solution of equation (2.11) for Uk 
with given Vk and prescribed /g), we also use linear finite elements for the 
discretization of Uk- In contrast to the data generation routine, we discretize 
the potential on a regular, structured grid which is usually much coarser than 
the grid used for the data generation. Later on, we shall assume that the phase- 
field Vk is also an element in the space of piecewise linear functions on the same 
underlying regular grid as for Uk ■ For the assembling of the stiffness matrix for 



( |2.11[ ), however, we replace the phase-field v^ by its L^-projection onto the space 
of functions which are piecewise constant on the triangles of the finite element 
space. A completely analogous procedure is applied for the solution of the adjoint 



equation (2.131 described in step 2 for the adjoint variable (pk- Note that both 



systems share the same stiffness matrix and that the right-hand side of (2.13) 
can be easily assembled using a slightly modified stiffness matrix. We shall use 
up to six different Cauchy data-sets for the reconstruction of the defect. The 
data-sets correspond to all possible combinations of pairs of electrodes where 
each electrode is located on a different side of the computational rectangle. We 
can use the same factorization of the stiffness matrix for all different right-hand 



sides of ( |2.11[ ) and ( |2.13 ). 

The calculation of the descent direction for the cost functional as described 
in step 3 requires another solution of an elliptic boundary value problem for the 
variable Svk — T^{DTe{vk))- As mentioned above, the update 6vk is discretized 
using linear triangular elements on a regular grid. To find dvk we have to solve an 
elliptic equations with system matrix defined by a discretization of the operator 
T : W* — ?> W. In our 2-dimensional test examples, we always set W — VF^'^(ri) 
and for any w* G W* we set T{w*) ~ v where v solves in a weak sense v — cAv — 
w* with some parameter c > and homogeneous Dirichlet boundary conditions. 
The choice of Dirichlet boundary conditions is motivated by the desire to keep 
the phase-field constantly at the value 1 on the boundary. The assembling of the 



right-hand side of the equation for 6vk is done by evaluating (2.14) for piecewise 
linear in all bases functions v. 

The projection required in step 4 is easily implemented for piecewise linear 
functions by thresholding the nodal values. Moreover, a suitable step-length for 
the update of the phase-field is found using an Armijo-type line search. We use 
a maximum number of five reduction steps for the correction of the step-length. 
Since each evaluation of the cost functional requires one solution of the state 
equation, we try to steer the step-size modification in a rather conservative way. 

Within this setup, the following numerical experiments have been performed. 
For all experiments, the phase-field parameter e was decreased in several steps 
from an initial value of e = 2 • 10"** down to e = 1 • 10~^ for the single-well 
potential and to £ = 2 • 10"^ for the double-well case. We run 2500 iterations of 
our algorithm in the single- well case and 1000 in the double- well case. Figure [l] 
shows the final phase-field together with the linear crack (as a white line) which 
was used for the data generation. We use all six available data-sets with electrode 
positions on (up/down), (left /right), (down/left), (up/left), (down/right), and 
(up/right) sides of the rectangle for the reconstruction and set the noise-level 
to zero. In this simple situation where the crack is located rather close to the 
boundary we obtain very good reconstruction of the crack location with the 
single-well approximation. 

In Figure [2] it is shown a comparison between reconstructions using 3 mea- 
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Figure 1: Reconstruction of a small linear crack with noise- free data. 

surements (left image) with electrode positions on (left /right), (left/up), and 
(right/up) edges and 6 measurements (right figure), again in the single- well 
case. It is notable that in the reconstruction with 3 data-sets the crack tips 
are accurately identified but the reconstructed crack is strongly curved which 
is probably due to the fact that we have no electrode located on the lower edge 
of the computational domain. In contrast the overall geometrical shape of the 
crack is reconstructed much better with 6 data-sets but the position of the crack 
tips is less accurate. In these two simulations we added one percent of normally 
distributed noise to Neumann and Dirichlet data. 




Figure 2: Comparison of reconstructions from 3 and 6 measurements. 

Figure [3] shows results for a situation with two cracks and different noise 
levels. Here we fixed the noise-level for the Neumann data to 1% for both exper- 
iments whereas the Dirichlet data were contaminated with 1% (left image) and 
5% (right image) of noise. We used three measurements (left /right), (left /up), 
(right/up) and the single-well potential. There is no big difference in the quality 
of the reconstructions. In both cases the placement of the smaller crack in the 
upper right corner is inaccurate and the larger crack in the lower left corner is 
curved. Nonetheless the convergence of the algorithm is not heavily effected by 
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the presence of (moderately strong) noise and the reconstructions are stable. 




Figure 3: Comparison with different noise levels. 

The next series of experiments presented in Figure HI shows the tendency 
of the single-well based algorithm to produce dendrite-like structures. In fact, 
the dendrite-shaped crack in the leftmost image is reconstructed quite well. The 
polygonal crack in the middle image is approximated by a cloth-hanger like 
structure which has a satisfactory data fit with a shorter overall length than the 
polygonal curve. Finally the cavity in the rightmost image is approximated by 
a one-dimensional structure which looks roughly like the skeleton of the cavity. 
In all these three experiments noise level is 1% for Neumann data and 5% for 
Dirichlet data and the three measurements (left/right), (left/up), (right/up) are 
used. 




Figure 4: Dendritc-like reconstructions with single-well potential. 

Figure [5] shows reconstructions obtained by using the double- well approxi- 
mation. As expected, the phase-field approximates the characteristic functions 
of one cavity (left image) and two cavities (right image) . In these two tests noise 
level is 1% for Neumann data. In the left image noise level for Dirichlet data is 
5% and the three measurements (left/right), (left/up), (right/up) are used. In 
the right image noise level for Dirichlet data is 1% but only one measurement, 
namely (left /right), is used. The overall location of the cavities is satisfactory, 
but the lower left quadrilateral is approximated by a non-convex shape. In this 
respect the experiment with the double-well potential for two cavities resem- 
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bles the results shown in Figure |2] where the lower left crack also has a strong 
tendency to bend inward. 

For our final numerical experiment, documented in Figure [SJ the double- 
well approach was used for the reconstruction of one-dimensional defects like 
the polygonal crack shown in the left image and the star-shaped crack shown 
on the right-hand side of the figure. In both cases the defect is approximated 
by a two dimensional structure. An interesting feature is the occurrence of a 
self-intersection of the boundary curve of the reconstructed defect in the case 
of the star-shaped crack. Also in these two final tests, noise level is 1% for Neu- 
mann data and 5% for Dirichlet data and the three measurements (left/right), 
(left/up), (right/up) are used. 




Figure 5: Reconstructions of cavities with double-well potential. 




Figure 6: Reconstructions of cracks with double-well potential. 

As a conclusion we can state that both algorithms give reconstructions of 
the defects with a satisfactory accuracy for an exponentially ill-posed problem. 
The algorithms show a quite stable behaviour in the presence of data noise. The 
single-well and double-well models develop the types of structures for which 
they are designed (one-dimensional for the single-well and two dimensional for 
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the double- well potential), so the single- well approach approximates cavities by 
dendrites and the double-well approach approximates cracks by cavities. The 
double-well approach looks more stable with respect to noise, is slightly less 
sensitive with respect to the adjustment of the phase-field parameter e and 
usually needs less iterations for convergence. This may be in accordance with 
the theory, in fact for the double- well case a convergence analysis is proved, 
whereas the single-well model we use is a modification of the one for which we 
have convergence results. Finally, it turned out to be important to update the 
phase-field parameter e adaptively during the algorithm. If the parameter s is 
chosen too small initially or decreased too fast, sharp interfaces develop too 
early, sometimes at incorrect locations, and the algorithm is not able to move 
well established interfaces to other locations. On the other hand, if the parameter 
e is decreased too much, the term containing the potential might prevail and 
not well established defects, usually the smaller ones, may disappear. 
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